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Abstract 

The Raise and Peel model is a recently proposed one-dimensional statistical model describing a 
fluctuating interface. The evolution of the model follows from the competition between adsorption 
and desorption processes. The model is non-local due to the possible occurrence of avalanches. At 
a special ratio of the adsorption-desorption rates the model is integrable and many rigorous results 
are known. Off the critical point, the phase diagram and scaling properties are not known. In this 
paper, we search for indications of phase transition studying the gap in the spectrum of the non- 
hermitian generator of the stochastic interface evolution. We present results for the gap obtained 
from exact diagonalization and from Monte Carlo estimates derived from temporal correlations of 
various observables. 
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I. INTRODUCTION 

The Raise and Peel model (RPM) |l| is a one-dimensional adsorption-desorption model of 
a fluctuating interface defined on a discrete lattice with restricted solid-on-solid rules. The 
adsorption process raises the interface while in the desorption process the top layer of the 
interface evaporates (it "peels" the interface) through non-local avalanches. The RPM shows 
self organized criticality 

as 

. The details of the dynamic will be presented in Subsection 

nn 

The continuum time limit of the model is characterized by a single coupling describing 
the ratio u of the interface adsorption/desorption rates. At u = 1, the model describes a 
c = logarithmic conformal field theory [4] with a dynamic critic exponent z — 1. It can 
also be analytically studied using its hidden Temperley-Lieb algebra [l[ [| . 

At this value of the adsorption/desorption ratio, many rigorous properties of the RPM 
are known about its equilibrium distribution, whose theoretical properties are related to 
alternating sign matrices and combinatorics p, Q . In the continuum limit in the temporal 
direction the ratio u can be any positive real number. In this paper we shall study the region 
< u < 1. 

Away from the Temperley-Lieb rates, much less is known about the RPM. In the prelim- 
inary analysis of Ref. the RPM at < u < 1 has been proposed to have a critical point 
u*, with < u* < 1, separating two phases. The phase at < u < u* is massive in the sense 
that the gap E\ of the non-hermitian generator of the stochastic temporal evolution tends 
to a finite positive value for large lattices. The phase in the region u* < u < 1 is instead a 
massless phase. The analysis of the behavior of Ei(u,L) as a function of u and the lattice 
size L has been crucial in making this conjecture. 

In this paper, we extend the analysis by pushing exact diagonalization up to the remark- 
able lattice length L = 26. We also propose a Monte Carlo estimate of E\ that we test up 
to L = 64. 

The analysis of the gap from these extended data provide information about the RPM 
phase diagram. Also, and most important, at the methodological level, it explores the finite 
size scaling properties of a non-local model discussing in a specific case the limits of the 
information that can be gained from this typical lattice sizes. 

The outline of the paper is the following. In Section II we recall the definition of the RPM 
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and discuss the numerical simulation technique in the (time) continuum limit. In Section 
III we recall the operatorial formulation of the model and summarize the conjectures about 
its phase diagram. Results on the exact diagonalization on larger lattices are showed in 
Section IV. Section V gives an estimate of the first excited level of energy E\ of this matrix 
while Section VI and VII present our Monte Carlo procedure to simulate the model and the 
analysis of the resulting data. Finally the Conclusions are summarized in Section VIII. 



II. THE RAISE AND PEEL MODEL 



A. Definition of the model 



The RPM belongs to the class of restricted solid-on-solid (RSOS) model. It is defined on 
a one dimensional lattice of size L + 1, with even L e 2N. The sites are numbered from 
to L, and at each site we assign a non negative height hi with RSOS constraints 

\h i+ x-hi\ = l, (2.1) 
h = h L = 0. (2.2) 

The space state Hl is discrete and its dimension is 

dimHL = (L/2 + lj(L/2r (2 - 3) 

To assign a stochastic discrete dynamics on the space state TCl, we first choose the 
two probabilities u a ,Ud G (0,1). In the language of Ref. [l|, they are associated to 
adsorption and desorption elementary processes. We then define the slope at site % as 
s% = h(hi+i — hi-i). The slope can only take values Sj G { — 1, 0, 1}. 

At each temporal step, we make the following stochastic move: 

1. choose at random a site 1 < % < L — 1 with uniform probability — 1); 

2. take the following actions (sub- moves) 
(a): if Si = and hi > hi-i, do nothing; 

(/3): if Si = and hi < hi-i, do hi — > hi + 2 with probability u a ; 
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(7): if Si = 1, with probability Ud, find j > i such that hj = hi and hk > hi for all 
i < k < j and make the transformation hk — > hk — 2 on all the internal sites 

% < k < j; 

(7'): if Si = —1, with probability Ud, find j < i such that foj = hi and /i fe > hi for 
all j < k < i and make the transformation hk — * /ifc — 2 on all the internal sites 
j < fc < i. 

It is easy to show that N 7 = Ny. 

In the following we shall be interested in the time continuum limit of this stochastic chain. 
This is discussed in the next paragraph where we also give details of how the continuous 
stochastic dynamics is practically simulated. 

B. (Time) continuum limit and numerical simulation of its stochastic dynamics 

The continuum limit in the temporal direction is taken at fixed L, by fixing the ratio 
u = u a /u d and sending u a ,Ud — > 0. The ratio u can be any positive real number. However, 
in this work, we shall consider mainly the region < u < 1, a restriction that we shall 
always assume if not said differently. In the conclusions, we shall briefly mention the case 
u > 1. 

In the limit Ud,u a — > 0, the lattice configuration is often frozen for many Monte Carlo 
iterations. Therefore, given a configuration s at time t, it would be very inefficient to perform 
the simulation for each time step t, t + 1, . . . ; it is better to compute the time interval T 
during which s stays unchanged, and the new configuration s' at time t + T. 

We group the sites in four classes according to the above ce, f3, 7 and 7' and denote by N a , 
N/3, etc., the numbers of sites of a certain class. We compute the total transition probability 
p = (u a Nj3 + 2udN 1 )/(L — 1). The probability distribution of T is proportional to (1 — p) T 
and can be sampled in the usual way by taking the logarithm of a uniform random number. 
Finally, we extract the site class f3, 7 or 7' with probability proportional to u a N/3, UdN 7 , 
uaNj, respectively, and pick the site % in the class with uniform probability. 
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III. OPERATORIAL DESCRIPTION OF THE RPM 



It is convenient to discuss the (continuum time limit of the) RPM and its dynamics 
with a Hilbert space formalism where we can easily rephrase all statements concerning the 
underlying stochastic process. We associate the vector \s) with each configuration s of the 
RPM. We denote by Ti the finite dimensional complex vector space having these states as 
an orthonormal basis 

(s\af) = 6 SS ,. (3.1) 
The subset TLs C H is the set of stochastic vectors of the form 

0<pi<l, = (3.2) 

i i 

They can be used to describe the state probability distribution at a certain time. 
We assign a special name to the state 

|a) = 5>>- (3.3) 

s 

It can be used to implement the sum over all states by taking scalar products. 

The evolution of the RPM is described by means of a non hermitian rate matrix H. The 
specific .natrix elects of Eb rL = 6 can be W^Q. 

An important property of H is the conservation of probability expressed by the relation 

(<t\H = 0. (3.4) 

Given H, we define the evolution operator U(t) = e~ tH . It has the following meaning 

U(t) s , s = (s'\U(t)\s) = Prob(s,to -> s',t + t). (3.5) 

In particular, if \ip) e Tis is a stochastic vector describing a statistical ensemble of states at 
time t Q , then the stochastic vector at time to + 1 is given by 

U(t)\1>)eH 8 . (3.6) 

The eigenvalues of H can be complex, and we sort them by increasing real part: 

H\<i/; n )=E n \ip n ), = E <ReE 1 <ReE 2 < (3.7) 



5 



It is convenient to introduce the expansion of the eigenvectors in terms of the single config- 
uration basis 

l^n) = J2^,s\s). (3.8) 

s 

The eigenvectors \tp n ) are not necessarily orthogonal. In the following we shall always assume 
that they are a basis of the state space. The condition KeEi > expresses the convergence 
property of the stochastic process described by H. Due to this condition we have 

lim U(t) s ' s = ^ov, (3.9) 

t^oo 

independently on the initial state s. In other words, every statistical ensemble converges 
asymptotically to the equilibrium distribution specified by the components of the lowest 
eigenvector. Due to the meaning of \tpo), it is convenient to choose the normalization of l^o) 
by requiring that it belongs to Hs 

J2^o,s = h or HVb) = l, (3.10) 

s 

in agreement with its probabilistic interpretation. 



A. The RPM away from the Temperley-Lieb point 

The special point u = 1 defines the so-called Temperley-Lieb (TL) rates. At this value of 
the adsorption/desorption ratio the RPM can be connected to the dense 0(1) or TL loop 
model 8| and many rigorous properties are known about its equilibrium distribution. 

In particular, at u — 1 the model is described by c = logarithmic conformal field theory 
(LCFT) with a dynamic critical exponent z = 1. 

Away from the Temperley-Lieb rates, much less is known about the RPM. The aim of 
this paper is precisely the investigation of the phase structure of the model with respect to 
the rate variable u. 

Following the analysis of Ref. the existence of possible critical points for u < 1 can 
be studied by computing the first excited level E\ as a function of the rate u and the lattice 



size L. 
of Ref. 



rom exact diagonalization up to L = 16 and finite size scaling analysis, the authors 



i] 



have claimed that the RPM exists in two phases when u < 1. The first phase is 
called massive and is realized for < u < u* . The second phase, for u* < u < 1 is a massless 
phase. 



As we mentioned in the Introduction, the naming refers to the different behavior of 
Ei(u,L) at large L. In a massive phase E% tends to a positive constant as L — > oo. In a 
massless phase, it tends to zero. At u — 1, we know rigorously that LE\ tends to a positive 
constant. 

r-i 

The finite size analysis of Ref. |l| studies the crossing of LEi(u,L) between the curves 
at the two values L and L + 2 with L < 16. Two crossings are observed. The smaller one 
appears to converge to some value u* around 0.5 as L increases. The larger crossing is less 
stable and it possible to conjecture that it tends to 1 at infinite L. Due to the small lattice 
sizes, these conclusions are quite preliminary and deserve a more detailed investigation. 

In the next Sections, we shall first discuss the outcome of an improved analysis based 
on exact diagonalization up to the size L = 26. Then, we shall illustrate a Monte Carlo 
method to compute E% and numerical data up to L = 64. These are less precise than exact 
diagonalization points, but can complement them, as we shall illustrate. 

To conclude this Section, we like to add a general comment. In principle, it is possible 
to study the critical behavior of the RPM by alternative methods, in particular by studying 
the size and u dependence of statistical averages over very large lattices. In our analysis we 
concentrated on the (difficult) measure of E\ with the precise aim of extending the above 
analysis. 



IV. EXACT DIAGONALIZATION ON LARGE LATTICES 

We have pushed the exact diagonalization of the RPM Hamiltonian up to the size L = 26, 
using the software package ARPACK 9]. A plot of LEi(u,L)) as a function of u for the 
various considered L is shown in Fig. (0). From that figure, we can see the existence of two 
crossings between the curves at L and L + 2. We shall denote them as u\ e f t (L) and it r ight(£) 
with Mi e ft(L) < M r ight {L) . We want to compute the extrapolations 

<rft, right = r lim Wtaft, right (£)■ (4.1) 

In Fig. we show the crossings u\ e ft(L), -u r i g ht(£) as functions of 1/L to test convergence. 
The left part of the figure shows u\ c f t {L). Data are accurately fit by the simple function 

0.6912 (1 - 8/L + 48/L 2 ). (4.2) 
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The integer coefficients are empirical, but the quality of the fit is remarkable with a x 2 — 
4.6 • 1CT 7 . 1 Our estimated extrapolation u* eit = 0.6912 suggests the existence of a critical 
point in infinite volume that separate different phases of the model. 

Of course, without further analytical insight, it is also perfectly possible that the above 
extrapolation is fake and that we are working on a lattice which is too small to extrapolate 
reliably to the true infinite size limit. We leave this kind of comments to the conclusions. 

The analysis is similar for « r i g ht(-^) shown in the right part of the figure. In this case, we 
also tried a polynomial in 1/L, but did not find a simple expression for the coefficients of the 
various powers of 1/L. The figure shows three curves obtained by fitting the first 8, 7, 6 data 
points with the highest L. As data points with smaller L are discarded, the extrapolated 
estimate for w* ight moves closer and closer to the special value 1. Thus our data support the 
conjecture u* ight = 1 with no additional non trivial critical point. The convergence in L is 
quite slow; this is not surprising, given the non-local nature of the dynamics. 

At size L = 26, the state space dimension is 742900. For computational reasons, it is not 
sensible to push the exact diagonalization beyond this point. Indeed, the convergence rate 
can be judged from the above figures and is not dramatically fast. 

For this reason, it is convenient to resort to Monte Carlo methods. In the next Sections, 
we shall first describe our theoretical approach, and then discuss the actual measurements 
of Ei on larger lattices. 

V. STOCHASTIC AVERAGES AND ESTIMATE OF Ei(u,L) 

The RPM is a Markov chain where the non hermitian Hamiltonian plays the role of a 
rate matrix. Our analysis of the phase structure of the RPM is based on the calculation of 
the first excited level Ei of this matrix. In this Section, we recall the basic formalism that 
allows to extract E\ from the time evolution of certain stochastic averages evaluated over 
the chain. 

1 Here and in the following, we denote by \ 2 the sum of the squares of the differences between the data 
points and the values of the fitting function. 
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A. Expansion of vectors 



We assumed that the eigenvectors of H are a basis of H. Hence, any state \(f>) can be 
expanded in terms of the eigenvectors of H 

|0) = ^c n |^„). (5.1) 

n>0 

The coefficient c can be determined. We apply U(t) to the above expansion and take the 
t — > +00 limit. Since E n > for n > 0, we find 

lim U(t)\<j>) = co\iPo). (5.2) 

Taking the scalar product with \a) and exploiting (a\U = (a\ as well as the chosen normal- 
ization of j^o) we conclude that Co = (cr\(f>), i.e., 

|0) = H0)^o) + ^c„|^). (5.3) 

n>l 

B. Time averages and their operator expression 

Let us consider various functions {Ai(s)}i<i<N of the RPM configurations. The stochastic 
A^-point averages times {U} starting from a given configuration s will be denoted by 

N 



(5.4) 

_i=i 

where the expectation is over all realizations of the stochastic process {s t }t>o with assigned 
initial condition. The equilibrium averages and correlations are independent on the initial 
state and are defined as 



E [A] = lim E any [A(t)\, (5.5) 
E[A(0)B(t)] = lim Eany [A(t)B(t + r)] , (5.6) 

(5.7) 



As a preparation to the calculation of stochastic time averages, we associate an operator 
to each function A(s) of the RPM configuration s. The rule is trivial. We associate to A(s) 
the operator defined by the following diagonal matrix elements in the \s) basis 

(s\A\s') = A(s)5 ss ,. (5.8) 



Consider now the stochastic average of a single A at time t, starting from the fixed 
configuration sq. By definition, we can write 

E S0 [A(s t )} = (a\AU(t)\s ). (5.9) 

From the above expansion we have 

|SO> = |^0> + 5> n |Vn> ( 5 - 10 ) 

n>l 

with certain coefficients c n . Thus, 

E S0 [A(s t )] = (a\A\i>o) + J2cne- tEn (a\A\^ n ). (5.11) 

n>l 

Of course, the asymptotic value is precisely the equilibrium average 

E[A] = (a\A\^ Q )=J2^ ,s. (5.12) 

s 

Another example is the self correlation of A with lag r. By definition, this is 

E [A(0)A(t)} = ^'Prob(s -> s')A s ^ s (5.13) 

ss' 

= Y, A s'U{r) s 'sA s ^ s = {a\AU{r)A\%) 

ss' 

Again, we have the expansion 

A\^) = b \^) + Y J K\^n) (5.14) 

n>l 

for certain coefficients b n . Then, the large r behavior of the connected correlation is given 
by 

E[A(0)A(r)]-E[A] 2 = £ c n (a\ A\^ )e~ E " T (5.15) 

n>l 

VI. MONTE CARLO SIMULATION AND ESTIMATE OF E x 

The Monte Carlo procedure consists of the following steps. 

A thermalized configuration s is generated by a straightforward implementation of the 
update procedure described in Subsection III Bl 

A set of b branches is generated by running b times the update procedure for a time 
duration r, each branch starting from s, with a different sequence of random numbers. 
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Along each branch, the values of the observables Ai(t) are averaged over time intervals of 
size At; the resulting values are used to estimate the expectation value of A,i{r At) with 
initial conditions s(0) = s: 

1 r(r+l/2)At 

A iir = A(r At) = — / (a\A U(t)\s). (6.1) 

At ./(r— 1/2) At 

We selected the following set of observables: 

L-l L-l-d 

A = h L/2 , A x =J2 hi > Ad + 2 = h t h i+d, < d < 14. (6.2) 

8=1 8=1 

We take as value and statistical error of Ai jt . the average and standard deviation over the 
branches at fixed time interval. We typically choose b = 10 8 . 

The resulting data for A iyr are fitted to the two-exponential formula 

A i)T = <$ + af exp(-/iP r At) + af exp(-/i^ } r At), (6.3) 

with fix < fi2, for r At > t m i n . The value and error of fix are estimated checking the stability 
of the fit vs. t m j n and the consistency of the best-fit values of fix between different observables. 

The procedure is repeated starting the branching process from several statistically- 
independent configurations s and the resulting values of fix are analyzed. If a configuration 
gives a value inconsistent with the others, it is removed. The value and error of fix are 
estimated from the remaining configuration by a standard statistical analysis; the error is 
multiplied by the scale factor s = a/ X 2 / n d.o.i if s > 1, i.e., if the values are not fully con- 
sistent within errors; we consider s 2 > 6 a warning sign of possible troubles. The resulting 
value and error of fix is directly an estimate of E\. 

A typical example is shown in Fig. El Some branch sets give extremely unstable results 
and are discarded from the analysis; it is likely that their initial configurations happen to 
have a very small superposition on |?/>o)- It turns out that A and Ax often give unstable fit 
results and are not consistent with each other and with other Af, therefore we exclude them 
from the analysis. 

VII. DATA ANALYSIS AND DISCUSSION 

We study the behavior of Ex vs. L at fixed u. We fit the exact data for 8 < L < 26 to 
one of the forms 

LEx(L) = c + ai L- u + a 2 L~ 2u (7.1) 
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(for the massless phase) and 

Ex(L) = c + ai exp(— pLiL) + a 2 exp(— ^L) (7.2) 

(for the massive phase). We discriminate the phase by the value of the \ 2 °f the fit and by 
consistency of the fit, extrapolated to large L, with Monte Carlo data. 

Two typical examples are shown in Figs.|3]and|Hlwhere we show the results of the analysis 
at the two points u = 0.4, and u = 0.9. The former value is clearly in the massive phase, 
the latter in the massless one. This is consistent with Fig. (J2J) at all considered u values, 
i.e. from 0.1 to 1.0 with 0.1 spacing. Of course, the discrimination in the critical region 
u ~ 0.7 is more difficult and less clear. In the massless phase, the determination of v is 
quite imprecise; it favors values v ~ 1/2, although, i.e., for u = 1 it is also compatible with 
v = 1. (In the massive phase, the "wrong" exponential fit often gives a negative v.) 

A combination of the exact diagonalization and Monte Carlo data is summarized in 
Fig. (jHJ). In each panel we have shown in log-log scale the ED data and the three points 
from Monte Carlo (L = 32, 48, and 64) versus L. 

The left panel contains data at u = 0.1 — 0.6. The ED data at 0.1-0.3 bend in the up 
direction suggesting convergence of E% to a finite positive limit as L —>■ oo. This conclusion 
is strongly supported by MC data. The bending is not visible at u = 0.4, 0.5, and 0.6, but 
again MC data suggest the same asymptotic behavior. 

The right panel of Fig. © shows the remaining values u = 0.7 — 1.0. In this case, the 
ED data bend in the down direction in agreement with an algebraic decay of Ei with the 
lattice size. This is confirmed by the MC points that extend smoothly this trend. In other 
words, at the explored lattice sizes, there are no signal for a plateau in E\. 

We also find remarkable that the special point u* ~ 0.69 singled out by the exact di- 
agonalization analysis remains a sort of separator point between two qualitatively different 
behaviors also after the inclusion of the Monte Carlo data points. 

VIII. SUMMARY AND DISCUSSION 

In summary, we have investigated the RPM away from its integrable point. First, we have 
extended the available results from exact diagonalization of the non local RPM Hamiltonian. 
As a further step, we have exploited Monte Carlo methods to obtain spectral information 
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on the Hamiltonian from the temporal decay of stochastic averages of various observables. 

As suggested in Ref. we have found evidence for a drastic change of behavior of 
the system at the intermediate adsorption/desorption rate ratio u* ~ 0.69. The analysis 
of data on lattices up to 64 sites long has confirmed rather clearly that the RPM is in a 
massive phase for u < u* where a plateau can be seen in the plot of E% vs. L. The region 
u* < u < 1, where our data show a steady decrease in E\, is more difficult to analyze. 
Monte Carlo simulations cannot exclude that Ei(u,L) reaches a positive asymptotic value 
at a typical lattice size growing rapidly as u — > 1. To exclude such a scenario it would be 
necessary to investigate the phase diagram of the RPM with other techniques able to reach 
very large lattices. We believe that a determination of E\ from Monte Carlo simulations 
is not feasible on lattices much larger than the length we explored. However, it would be 
possible to measure static averages of typical non local quantities, like cluster distributions 
in the interface configurations, with good accuracy on quite larger lattices. Preliminary 
results along this line seem to suggest the existence of a single massive phase for all u < 1 2 . 
If such a scenario were confirmed, the Raise and Peel model would be an interesting 
example of very slow convergence to criticality in the infinite volume limit, presumably 
related to the non-locality of the model. We emphasize again that the RPM is non-local 
at generic u because of the occurrence of avalanches. This means that the Hamiltonian H 
does not describe a lattice model with local interaction and it is not straightforward to 
apply field theory techniques. Only at the TL point u — 1, the hidden loop algebra permits 
a straightforward and clean analytical investigation. Off integrability, some better degrees 
of freedom must be identified for an effective description of the avalanche dynamics. 



Apart from these speculations, we believe that our proposal for a direct measure of E\ 
is valuable since the gap is often an important quantity to be measured on critical discrete 
models. In principle indeed, the strategy described in this paper could be applied to more 
involved versions of the RPM. In particular, it would be interesting to analyze the role of 
boundary conditions as discussed in Ref. [10]. 

2 V. Rittenberg, private communication 
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FIG. 1: Exact diagonalization data for the function LE\(u, L). The various lines shows the results 
at sizes L = 8, 10, . . . , 24, 26. 
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FIG. 5: Ei vs. L for u = 0.9. Both x 2 and MC data clearly favor the massless phase. 



18 



e 

II 

o 

H 

I 

H 

o 

B' 



CO 
CO 

> 



i — i 
P 



X 

o 
^+ 

o 

TO 
O 
P 

S' 

& 

o 

g 
a 

o 



O 
cT 



o 
i-i 



CO 



B' 

i 

cT 

Crq 

CO 

o 

CD 



0.4 

0.2 

0.0 E- 
-0.2 
-0.4 
-0.6 
m -0.8 E- 



-1.0 
-1.2 
-1.4 
-1.6 E- 
-1.8 
-2.0 



i ii in 



o. 



J °o c 





o 


■o u = 0.1 




□ ■ 


□ u = 0.2 




o- 


■O u = 0.3 




A 


A u = 0.4 




< ■ 


■■< u = 0.5 




V ■ 


■v u = 0.6 



i4 



I I I I I I I I I I I I I I I I I I I 



2.5 



3.5 
logL 



4.5 



0.5 



0.0 



-0.5 



-1.0 



-1.5 



-2.0 



4 l l l I l l l l | l l l l 

_S>. 1 
PY-. 
-QV\A 

•• o . 

\b\A 

OA*. A 

'•P.Va 



I I I I I I I I I 



■ 


o 


u = 0.7 


□■■ 


■■□ 


u = 0.8 


o- 


o 


u = 0.9 


A ■ 


A 


u= 1.0 



Qb<> 

oh 
b 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



2.5 



3.5 
logL 



4.5 5 



